/*************************************************************************
Purpose: Loads map output from GAIA and modifies it in some manner. 
         Usually used with unstructured maps, but latlon maps are accepted
         by some choices.
**************************************************************************
Copyright 2010 Bryan Killett

This file is part of GAIA.

GAIA is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

GAIA is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with GAIA. If not, see <http://www.gnu.org/licenses/>.
*************************************************************************/
//Header file "definitions.hpp" declares many functions, constants, datatypes.
//In order to use the functions declared in definitions.hpp, must link
//with object code that includes functions.cpp.
#include "definitions.hpp"

int main(int argc, char *argv[]){  
  //Variable declarations- every variable type that ends in "_s"
  //is defined in GAIAFOLDER/c/definitions.hpp.
  long i,j;//Counters.
  string input_file(gaiafolder()); //Filename of options file created in IDL.
  input_s input;
  ancillary_s ancillary;
  grid_s region;
  grid_s grid, grid2;
  model_s model;
  results_s results1,results2,final_results;
  plot_options_s plot_options,plot_options2,temp_options;
  insitu_s insitu;//For tide gauge data, etc.
      
  string base;//This string forms the base part of the output filename.

  string input_folder1,input_folder2;
  
  //////////////////////////////////////////////////////////////////////////
  //End of variable declarations, start of initialization/loading section.
  //////////////////////////////////////////////////////////////////////////
  
  //Each program should manually set this value.
  //See definitions.hpp for more details.
  results1.options.output_choice = 1;
  results2.options.output_choice = results1.options.output_choice;
  final_results.options.output_choice = results1.options.output_choice;

  //Controls how old output files will be modified to create new output files.
  //Old output files are found in input_folder1,2.
  //Both folders need to have the same number of files and grid points and
  //locations of grid points.
  //Case 1 doesn't require its folder to have output_filenames.txt, but
  //other cases do.
  //New output files will be saved to folder specified in argv options.
  //1   - Compile many files like region_00001_plot_data_0001.txt together.
  //      No output_filenames.txt needed- folder1 is searched for files.
  //11  - Subtract the files in the 2nd listed output folder from the 1st.
  //      Don't do any sin/cos conversion- just subtract the maps.
  //12  - Add the files in the 1st and 2nd listed output folders, again very crudely.
  //21  - Subtract the files in the 2nd listed output folder from the 1st.
  //      Do sin/cos conversion and convert back to amp/phase if both folders
  //      are in that format.
  //22  - Add the files in the 1st and 2nd listed output folders in sin/cos format.
  //      Do sin/cos conversion and convert back to amp/phase if both folders
  //      are in that format.
  //23  - Result = files in the 1st listed output folder divided by those in the 2nd
  //      Do sin/cos conversion and convert back to amp/phase if both folders
  //      are in that format.
  //31  - Convert files in 1st folder from amp/phase to sin/cos.
  //32  - Convert files in 1st folder from sin/cos to amp/phase.
  //41  - Multiply output in files in 1st folder by parameter1.
  //101 - Convert from J2000 epoch to Greenwich phase lags.
  //102 - Convert from Greenwich phase lags to J2000 epoch.
  //103 - Convert from input.epoch to J2000 epoch
  //104 - Convert from J2000 epoch to input.epoch
  //201 - Calculate avg land/ocean sin/cos coefs and std devs.
  //301 - Convert from A*cos(wt) + B*sin(wt) to a*cos(wt+lon) + b*cos(wt-lon)
  //302 - Convert from A*cos(wt) + B*sin(wt) to a*sin(wt+lon) + b*sin(wt-lon)
  int modification_choice = 201;
  //load_output wants the folder names relative to GAIAFOLDER.
  //input_folder1 = "output/temp/";
  input_folder1 = "fortran/paper/Figs/real/20110324-012805 - 7 years fixed orbit verif last minus verif first 4M, 230km gr supp3, q1,n2,mf,mm/";
  input_folder2 = "output/maps/20110324-013540 - FES plus inversion, q1,n2,mf,mm/";
  
  //This parameter is used in several modification cases.
  double parameter1=pi;

  //Load options from file in folder 1.
  input_file.append(input_folder1); input_file.append("options.txt");
  load_options(input_file, input, region.options, ancillary, grid, model.options, results1, temp_options, insitu, 1);

  //Use current time and command line args to create output folder.
  create_output_folder(argc,argv,plot_options);
  
  //Sometimes I need to store the outputfolder in another variable temporarily
  string output_folder(plot_options.outputfolder.c_str());

  cout<<"Loading output1 from folder1 "<<gaiafolder()<<input_folder1<<endl;
  if(modification_choice != 1) cout<<"Loading output2 from folder2 "<<gaiafolder()<<input_folder2<<endl;
  cout<<"Saving modified output to "<<plot_options.outputfolder<<endl;

cout<<"!!!NEED TO CP OPTIONS FILE AND MISC_results AND OTHER FILES TO NEW DIR!"<<endl;
  
  //Load ancillary data (topography, temperatures, etc) from disk.
  load_ancillary(ancillary, 1);

  //load_output() expects a grid to already be defined.
  define_grid(grid, ancillary, results1);

  //I preset the input folder to 1 because that's what most cases do anyway.
  plot_options.inputfolder = input_folder1;

  switch(modification_choice){
    ///////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////
    default:cout<<"WARNING!!! Didn't recognize modification_choice"<<modification_choice<<endl;
    ///////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////
    case 1:{//Compile partial output from parallel region run.
      region = define_region(region.options, grid, 1);
      //Need options because h_size comes up a lot in parallel loading.
      final_results.options = results1.options;
      //find_and_load needs this variable defined correctly.
      plot_options.outputfolder=gaiafolder();//output folder needs to be absolute
      plot_options.outputfolder.append(plot_options.inputfolder);
      find_and_load_parallel_region_output(region, plot_options, final_results, grid);
      //Now return outputfolder to its correct value for writing.
      plot_options.outputfolder=output_folder;
      break;
    }
    ///////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////
    case 11:{//Subtract two output folders crudely.
      //First load both folders into their respective results structures.
      cout<<"Loading the first folder."<<endl;
      load_output(results1,plot_options,grid);
      plot_options2.inputfolder = input_folder2;
      cout<<"Loading the second folder."<<endl;
      load_output(results2,plot_options2,grid);
      //Quick check- do both have same number of parameters and grid pts?
      if(results1.outputs.size() != results2.outputs.size()) cout<<"First folder has "<<results1.outputs.size()<<" params while second folder has "<<results2.outputs.size()<<" params."<<endl;
      if(results1.outputs.at(0).size() != results2.outputs.at(0).size()) cout<<"First folder has "<<results1.outputs.at(0).size()<<" grid pts while second folder has "<<results2.outputs.at(0).size()<<" grid pts."<<endl;
      
      //Copy titles and units from 1st folder.
      final_results.titles = results1.titles;
      final_results.units = results1.units;
      
      //Add Diff'd to all titles.
      for(i=0;i<(long)final_results.titles.size();i++)
        final_results.titles.at(i) = "Diff'd "+final_results.titles.at(i);
      
      //Resize final_results using 1st folder.
      final_results.outputs.resize(results1.outputs.size());

      //The ranges of the second indices are all equal to the # gr pts in rg
      for(i=0;i<(long)final_results.outputs.size();i++)
        final_results.outputs.at(i).resize(results1.outputs.at(0).size(),0.0);
      
      //Now compute 1 - 2, store in final.
      //Loop through parameters.
      for(i=0;i<(long)results1.outputs.size();i++)
        //Loop through grid points to copy output, but treat phases different!
        if(results1.options.output_type.at(i)==5){//Phase detected.
          for(j=0;j<(long)results1.outputs.at(0).size();j++)
            final_results.outputs.at(i).at(j) = restrict_degrees(results1.outputs.at(i).at(j) - results2.outputs.at(i).at(j));
        }
        else{
          for(j=0;j<(long)results1.outputs.at(0).size();j++)
            final_results.outputs.at(i).at(j) = results1.outputs.at(i).at(j) - results2.outputs.at(i).at(j);
        }  

      break;
    }
    ///////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////
    case 12:{//Add two output folders crudely.
      //First load both folders into their respective results structures.
      cout<<"Loading the first folder."<<endl;
      load_output(results1,plot_options,grid);
      plot_options2.inputfolder = input_folder2;
      cout<<"Loading the second folder."<<endl;
      load_output(results2,plot_options2,grid);
      //Quick check- do both have same number of parameters and grid pts?
      if(results1.outputs.size() != results2.outputs.size()) cout<<"First folder has "<<results1.outputs.size()<<" params while second folder has "<<results2.outputs.size()<<" params."<<endl;
      if(results1.outputs.at(0).size() != results2.outputs.at(0).size()) cout<<"First folder has "<<results1.outputs.at(0).size()<<" grid pts while second folder has "<<results2.outputs.at(0).size()<<" grid pts."<<endl;
      
      //Copy titles and units from 1st folder.
      final_results.titles = results1.titles;
      final_results.units = results1.units;
      
      //Add Summed to all titles.
      for(i=0;i<(long)final_results.titles.size();i++)
        final_results.titles.at(i) = "Summed "+final_results.titles.at(i);
      
      //Resize final_results using 1st folder.
      final_results.outputs.resize(results1.outputs.size());

      //The ranges of the second indices are all equal to the # gr pts in rg
      for(i=0;i<(long)final_results.outputs.size();i++)
        final_results.outputs.at(i).resize(results1.outputs.at(0).size(),0.0);
      
      //Now compute 1 + 2, store in final.
      //Loop through parameters.
      for(i=0;i<(long)results1.outputs.size();i++)
        //Loop through grid points to copy output.
        for(j=0;j<(long)results1.outputs.at(0).size();j++)
          final_results.outputs.at(i).at(j) = results1.outputs.at(i).at(j) + results2.outputs.at(i).at(j);
 
      break;
    }
    ///////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////
    case 21: //Subtract   two output folders in sin/cos format.
    case 22: //Add        two output folders in sin/cos format.
    case 23:{//Divide     two output folders in sin/cos format.
      //First load both folders into their respective results structures.
      cout<<"Loading the first folder."<<endl;
      load_output(results1,plot_options,grid);
      plot_options2.inputfolder = input_folder2;
      cout<<"Loading the second folder."<<endl;
      load_output(results2,plot_options2,grid);
      //Quick check- do both have same number of parameters and grid pts?
      if(results1.outputs.size() != results2.outputs.size()) cout<<"First folder has "<<results1.outputs.size()<<" params while second folder has "<<results2.outputs.size()<<" params."<<endl;
      if(results1.outputs.at(0).size() != results2.outputs.at(0).size()) cout<<"First folder has "<<results1.outputs.at(0).size()<<" grid pts while second folder has "<<results2.outputs.at(0).size()<<" grid pts."<<endl;

      convert_results_from_ampphase_to_sincos(results1);
      convert_results_from_ampphase_to_sincos(results2);
      
      if(modification_choice == 21){
        cout<<"Calculating folder1 - folder 2"<<endl;
      }
      else if(modification_choice == 22){
        cout<<"Calculating folder1 + folder 2"<<endl;
      }
      else if(modification_choice == 23){
        cout<<"Calculating folder1 / folder 2"<<endl;
      }

      //Copy everything from 1st folder.
      final_results = results1;
      
      //Loop through parameters.
      for(i=0;i<(long)results1.outputs.size();i++){
        //Loop through grid points.
        for(j=0;j<(long)results1.outputs.at(0).size();j++){
          if(modification_choice == 21){
            final_results.outputs.at(i).at(j) = results1.outputs.at(i).at(j) - results2.outputs.at(i).at(j);
          }
          else if(modification_choice == 22){
            final_results.outputs.at(i).at(j) = results1.outputs.at(i).at(j) + results2.outputs.at(i).at(j);
          }
          else if(modification_choice == 23){ 
            if(fabs(results2.outputs.at(i).at(j)) > 0.1)
              final_results.outputs.at(i).at(j) = results1.outputs.at(i).at(j) / results2.outputs.at(i).at(j);
            else
              final_results.outputs.at(i).at(j) = 0.0;
            
            //Division is intended to look at residual GRACE divided by FES, so discard > 100%.
            if(fabs(final_results.outputs.at(i).at(j)) > 1.0) 
              final_results.outputs.at(i).at(j) = 0.0;

          }
          else cout<<"!!!!WARNING!!! ERROR DETECTED!!!"<<endl;
        }
      }

      //If both were originally amp/phase, convert sin/cos output back to amp/phase
      if(results1.conversion_to_sincos == 2 and results2.conversion_to_sincos == 2) convert_results_from_sincos_to_ampphase(final_results);
      
      //Now modify all titles so they won't be over-written by sin/cos conversion.
      for(i=0;i<(long)final_results.titles.size();i++){
        if(modification_choice == 21){
          final_results.titles.at(i) = "Diff'd "+final_results.titles.at(i);
        }
        else if(modification_choice == 22){
          final_results.titles.at(i) = "Summed "+final_results.titles.at(i);
        }
        else if(modification_choice == 23){
          final_results.titles.at(i) = "Div'd "+final_results.titles.at(i);
          //Also, division changes units.
          final_results.units.at(i) = "#";
        }
      }
      break;
    }
    ///////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////
    case 31:{//Convert files in 1st folder from amp/phase to sin/cos.
      cout<<"Loading the folder."<<endl;
      load_output(final_results,plot_options,grid);
      convert_results_from_ampphase_to_sincos(final_results);
      break;
    }
    ///////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////
    case 32:{//Convert files in 1st folder from sin/cos to amp/phase.
      cout<<"Loading the folder."<<endl;
      load_output(final_results,plot_options,grid);
      convert_results_from_sincos_to_ampphase(final_results);
      break;
    }
    ///////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////
    case 41:{//Multiply output in files in 1st folder by parameter1.
      //First load folder.
      load_output(final_results,plot_options,grid);
      
      //Now multiply by parameter1.  Loop through parameters:
      for(i=0;i<(long)final_results.outputs.size();i++)
        //If this isn't a phase, Loop through grid points to copy output.
        if(final_results.options.output_type.at(i) != 5) for(j=0;j<(long)final_results.outputs.at(0).size();j++)
          final_results.outputs.at(i).at(j) *= parameter1;
 
      break;
    }
    ///////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////
    case 101:
    case 102:
    case 103:
    case 104:{
      cout<<"Loading the folder."<<endl;
      load_output(final_results,plot_options,grid);
      region = define_region(region.options, grid, 1);
      load_input(input,grid,region);
      final_results.options.adjust_phases_choice = modification_choice-100;
      adjust_phases(input,final_results);
      break;
    }
    ///////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////
    case 201:{//Calculate avg land/ocean sin/cos coefs and std devs.
      cout<<"Loading the folder."<<endl;
      load_output(final_results,plot_options,grid);
      estimate_land_tide(final_results,grid);
      break;
    }
    ///////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////
    case 301: //Convert from A*cos(wt) + B*sin(wt) to a*cos(wt+lon) + b*cos(wt-lon)
    case 302:{//Convert from A*cos(wt) + B*sin(wt) to a*sin(wt+lon) + b*sin(wt-lon)
      cout<<"Loading the folder."<<endl;
      load_output(results1,plot_options,grid);
      convert_results_from_ampphase_to_sincos(results1);
      
      //Copy everything from 1st folder.
      final_results = results1;
      
      //Change all titles.
      for(i=0;i<(long)final_results.titles.size();i++)
      
      //Loop thru outputtypes (periods count twice)- only affect amp/phases.
      for(j=0;j<(long long)final_results.options.output_type.size();j++){
        if(final_results.options.output_type.at(j)==2){//Cosine coef.
          if(modification_choice == 301){
            final_results.titles.at(j)   = "Cos +lon "+final_results.titles.at(j);
            final_results.titles.at(j+1) = "Cos -lon "+final_results.titles.at(j+1);
          }
          else if(modification_choice == 302){
            final_results.titles.at(j)   = "Sin +lon "+final_results.titles.at(j);
            final_results.titles.at(j+1) = "Sin -lon "+final_results.titles.at(j+1);
          }
          
          //Loop through grid points.
          for(i=0;i<(long)final_results.outputs.at(j).size();i++){
            double cos_coef = final_results.outputs.at(j).at(i);
            double sin_coef = final_results.outputs.at(j+1).at(i);
            double lon = deg2rad*grid.lon.at(i);
            if(modification_choice == 301){
              //Replace "A" with "a".
              final_results.outputs.at(j).at(i)   = 0.5*(cos_coef/cos(lon) - sin_coef/sin(lon));
              //Replace "B" with "b".
              final_results.outputs.at(j+1).at(i) = 0.5*(cos_coef/cos(lon) + sin_coef/sin(lon));
            }
            else if(modification_choice == 302){
              //Replace "A" with "a".
              final_results.outputs.at(j).at(i)   = 0.5*(cos_coef/sin(lon) + sin_coef/cos(lon));
              //Replace "B" with "b".
              final_results.outputs.at(j+1).at(i) = 0.5*(cos_coef/sin(lon) - sin_coef/cos(lon));
            }
          }
        }
      }//End of loop through output_type.
      break;
    }
    ///////////////////////////////////////////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////////////////////////////////////////
  }//End o switch-case-break.
  
  //Set the output filename base so that these output files are final.
  base = "plot_data";
  
  //Now that the analysis has been performed, it's time to output the data
  //in a format that the IDL program grace_input_map.pro can turn into a 
  //plot.  Last variable is 1 because this is the final write_output.
  write_output(base, final_results, grid, plot_options,1);

  return 0;
}
